Week 12 - Lecture Demo: Spectral Indices, Landsat ARD, Interpolation

Landsat ARD Coverage Extent

Concepts & Themes:

This week’s lab will feature exploration of the following 3 themes:

  1. Application of Spectral Indices to Landsat ARD imagery
  2. Utilization of Surface Temperature produce from Landsat ARD
  3. Raster Interpolation

Lecture Demo Data:

Part I Lab Steps - Spectral Indices

What is a spectral index?

  • A spectral index is a mathematical equation that is applied on the various spectral bands of an image per pixel.

  • The most common mathematical formulas that are used is the normalized difference:

(Bx - By)/(Bx + By)

  • In practical terms, this is the difference between two selected bands normalized by their sum. This type of calculus is very useful to minimize (as much as possible) the effects of illumination (shadows in mountainous regions, cloud shadows, etc) and enhance spectral features that are not visible initially. In other words, the calculation ‘normalizes’ the image, resulting in scaled values and a easily understood color ramp.

  • To Start, the input data for the lecture demo will be the same as that of the assignment itself, accessed through Earth Explorer:

Landsat C2 US ARD

  • Utilize the horizontal and vertical position of the ARD Grid to produce results just for that one grid space. For the lecture demo, the urban area of interest is Springfield, Missouri which has a Horizontal ID of 18, and Vertical ID of 11:

Horizontal/Vertical Position ID of AOI ARD Tile

  • Next, the following three conditions are met from the returned results:
  1. Cloud-free over Area of Interest - AOI - Springfield, Missouri
  2. The tile preview shows full pixel coverage over AOI
  3. The ARD tile date is within ‘the growing season’, i.e. approximately late May through early August. This is necessary as the analysis will feature spectral indices for vegetative cover, and we need the coverage to occur during the general growing season, not dormant seasons.

This Tile is cloud-free, full coverage over the AOI, and is within the growing season - 6/20/21

  • The Download features two products - the Surface Reflectance Values; and secondly the Surface Temperature values:

Surface Reflectance Bundle

Surface Temperature Bundle

  • To Start, the Surface Reflectance Values band 1 - 7 have been composited into SR_springfield_composite.tif located in the Part I directory. We can set the band combination to 5/4/3 to see the vegetative cover in/around the city relative to impervious surfaces:

Preview Imagery in Composite 5/4/3

  • Unlike the previous assignment whereby we used a composite image as the base for image classification, to create a spectral index, we need to access the original bands in the imagery. To start, we will import the bands necessary for the first spectral index of the lecture demo - NDVI. To follow, the bands shown in the layers panel along with the NDVI formula:

NDVI = (NIR Band - Red Band / (NIR Band + Red Band)

  • As translated to the Landsat 8 band numbers:

NDVI = (Band 5 – Band 4) / (Band 5 + Band 4)

Utilize Either the Composite Bands 4 and 5, or Import Them Separately

  • The following formula utilizes the bands as separate of the composite imagery SR_springfield_composite.tif:

("LC08_CU_018011_20210620_20210703_02_SR_B5@1"-"LC08_CU_018011_20210620_20210703_02_SR_B4@1")/("LC08_CU_018011_20210620_20210703_02_SR_B5@1"+"LC08_CU_018011_20210620_20210703_02_SR_B4@1")

  • Review the results with symbolization applied for values within the NDVI Index:

NDVI low > high Symbolization

  • The following graphic is useful to interpret the resulting NDVI values relative to vegetative cover:

NDVI Values Interpretation

  • Next, the same process of band combinations will be made, but this time for the The Normalized Difference Built-up Index (NDBI), which uses the the NIR and SWIR bands to emphasize manufactured built-up areas:

NDBI = (SWIR - NIR) / (SWIR + NIR)

  • As translated to the Landsat 8 band numbers:

NDBI = (Band 6 – Band 5) / (Band 6 + Band 5)

  • The following formula utilizes the bands as separate of the composite imagery SR_springfield_composite.tif:

("LC08_CU_018011_20210620_20210703_02_SR_B6@1" - "LC08_CU_018011_20210620_20210703_02_SR_B5@1") / ("LC08_CU_018011_20210620_20210703_02_SR_B6@1" + "LC08_CU_018011_20210620_20210703_02_SR_B5@1")

  • Review the results with symbolization applied for values within the NBDI Index. Notice the very close coorespondence with high built-up values with low NDVI values:

NBDI vs NDVI

Part II Lab Steps - Land Surface Temperature

  • Working with remote sensed imagery in a desktop environment typically involves much pre-processing work devoted to conversions and corrections. This paper does an excellent job of summarizing these pre-processing tasks for Landsat imagery:

  • A survival guide to Landsat preprocessing

  • One of the more onerous analysis tasks is the pre-processing of imagery for Land Surface Temperature - a topic that is trending for all the wrong reasons due to climate change coupled with social, demographic and built environmental inequities, especially in large cities.

  • Utilizing the ARD dataset for Landsat imagery, a good portion of pre-processing work has been completed by Landsat/USGS.

  • Currently in tar.gz compression, the directory is unzipped, and the ST raster is isolated. This raster is imported into QGIS via the Data Source Manager as a typical raster product:

The Surface Temperature ST Band

Surface Temperature (ST) – Represents the temperature of the Earth’s surface in Kelvin (K).

  • Also noted in the guide is a scale factor that is utilized for the product.

Kelvin Values

  • ST Band 10 in QGIS, Layers Panel:

Note the Kelvin Values 28804 - 49411

  • Here the digital numbers need to be scaled accordingly. If the value is to be expressed in Celcius or Fahrenheit, a further conversion needs to take place via the raster calculator. This will be accomplished in two steps:

    • Scale the raster values in raster calculator

    • Apply the following temperature conversion where K is the scaled raster values and F is the Fahrenheit temperature measurement:

    • F = 1.8(K - 273) + 32

  • To Start, open the Raster Calculator and apply the following scale factor and offset value as noted in the product metadata, exporting the product to the data folder as st_scaled:

    • "LC08_CU_018011_20210620_20210703_02_ST_B10@1" * 0.00341802 + 149

Scaled + Offset Kelvin Values

  • Next, apply the temperature conversion based on the formula for kelvin to fahrenheit units, and output as st_fahrenheit

    • 1.8 * ("st_scaled@1" - 273) + 32

  • Review the resulting raster values, Kelvin vs Fahrenheit:

Kelvin vs Fahrenheit

  • Utilize QuickMapServices plugin to compare high and low ST values relative to landcover:

ST values relative to Landcover

  • Finally, compare high resolution base imagery to NDVI to ST values:

Comparison of landcover to NDVI to ST values

Part III Lab Steps - Interpolation & Density Surfaces:

  • The data input for Part II has been derived from the Vehicle Crash feed at Open Data NY. While a bit gruesome, this is a good technical input as there is no discernable pattern to the data on load to QGIS (316623 crashes as of 11/8/2019). A density surface can uncover the ‘pattern to the data’ beyond the points features themselves:

Part III Data Structure

Part III Input Data Loaded to QGIS Layers Panel

  • Its important to note at this juncture that there are two primary reasons for surface interpolation. First, like the current mapping, a pattern to the data beyond the feature input is sought. Second, and most related to climatology datasets used this week in the assignment, is the intent to ‘fill’ in missing data and ‘smooth’ data across geographic expanse. This is the situation where values between location A and location B need to be ‘filled’ in. In effect, surface interpolation gives up the specificity of input location precision to get back better predictive values across spaces - including those ‘between’ known data input locations and their values.

  • Prior to running a density surface, a defined search radius needs to be determined. This is the euclidean distance from a individual cell in the resulting surface outwards to defined a ‘search distance’. This distance can be derived by a formula or it can be ‘eyeballed’. Based on a formula calculation, the optimal distance for the dataset was found to be 582 feet.

  • To learn more about how to derive a optimal distance - known as hopt - the following video outlines steps to do so.

    • Video Tutorial for Optimal Bandwidth (hopt) calculation for a radius distance for kernel density estimation: Determine Optimal Bandwidth (hopt)

  • To Start the process, there are two build-in tools and also a plug-in tool. Here we will use the Heatmap (Kernal Density Estimation) tool:

Heatmap (Kernal Density Estimation) Tool

Input Parameters

  • Input is the points feature MVC_11_8_2019

  • Radius = 582

  • Pixel Size X and Y = 100 (note: the feature is projected as NYSP, units Feet; thus all measurments are in square feet or linear feet, depending on the measurement type.

  • Output = Raw - this results in a estimated count per pixel of accidents

  • Kernel Shape = Quartic - default

  • There will be no weights; that is, just the locations of the accidents, not their variables, will be estimated.

  • View result:

Interpolation Result

  • Next, symbolize the results based on a spectral color ramp (invert to set high values to red, not blue). Further, cut the statistical range by setting the Cumulative cut count to 2 - 98 % to cut out any values below zero:

Symbolization

In Review of the resulting surface, apply 50% transparency and drape the kernel raster over a satellite base map. Here the ESRI Satellite base is utilized:

Review - Small Scale

  • In Review of the various hotspots throughout the kernel density result, its now clear that density of accidents is spatially consistent with bridges and tunnels throughout Manhattan. Without the density surface, this concentration pattern is largely missing from the input points feature itself. Thus we get a quantitative surface with cell values of accidents and and insightful visualization:

  • Lower Manhattan:

Review - Large Scale

  • Midtown:

Review - Large Scale

  • Uptown near GWB:

Review - Large Scale